//      From Unemployment to Self-employment: An Evaluation of Self-Employment Assistance Programs
/*      2023, Sumudu Kankanamge, Alexandre Gaillard
 
 
 
 /********************************** MAIN **************************************/
#include "parameter.h"
#include "tools/useful.cpp"
#include "tools/tauchen.cpp"

// Optimisation //
#include "optimization_various/SMMfun.cpp"                  /// to compute SMM distance///
#include "optimization_various/readinput.cpp"               /// to compute SMM distance///]
#include "optimization_various/zbrentNEW.cpp"               // to find minimum sequentially ///
#include "optimization_various/mymngolden.cpp"              // to find minimum /////
#include "optimization_various/myCRSwithoutparams.cpp"      // To find best params ///


// function entrepreneur's production and profit
#if OPT_labor_demand == 0
double labor_demand(double zz, int tgridfoc, double kk){return(0.0);}
double prodE(double zz, int tgridfoc, double kk){return((zz)*(gmapping[tgridfoc])*pow((kk),nupar));}
double profitE(double zz, int tgridfoc, double kk){return((((prodE((zz),tgridfoc,(kk)) - deltapar*(kk)))));}
#endif
#if OPT_labor_demand == 1
double labor_demand(double zz, int tgridfoc, double kk){
    double labor_tmp = pow((wstar/(zz*gmapping[tgridfoc]*nupar*(1.0-gammapar)))*(1/pow((kk),nupar*gammapar)),(1/((1.0-gammapar)*nupar-1.0)))-labor_supply_ent[tgridfoc];
    return(max(0,labor_tmp));
}
double prodE(double zz, int tgridfoc, double kk){
    return((zz)*(gmapping[tgridfoc])*pow((kk),nupar*gammapar)*pow((labor_supply_ent[tgridfoc]+labor_demand(zz,tgridfoc,kk)),(1.0-gammapar)*nupar));
}
double profitE(double zz, int tgridfoc, double kk){
    return((((prodE((zz),tgridfoc,(kk)) - labor_demand(zz,tgridfoc,kk)*wstar - deltapar*(kk)))));
}
#endif


// Self-employment Assistance Program.
double bstarfun(double returns, int tgridfoc) {
    double bstar = 0.0, btheta = 0.0;
    
#if OPT_policy_exp == 2
    btheta = min((1.0-taxrate)*rhostar*wstar*prod[0]*parascale_DRI,(1.0-taxrate)*rhostar*wstar*prod[tgridfoc]);
#endif
#if OPT_policy_exp != 2
    btheta = (1.0-taxrate)*rhostar*wstar*prod[tgridfoc];
#endif

    if(returns <= 0.0) {bstar = btheta + returns;}
    if(returns > 0.0 & returns < (btheta)/DRIpar){bstar = btheta + (1 - DRIpar)*returns;}
    if(returns >= (btheta)/DRIpar){bstar = returns;}

    
#if OPT_policy_exp == 3
    bstar = (1.0-taxrate)*rhostar*wstar*prod[tgridfoc] + returns;
#endif
#if OPT_policy_exp == 5
    bstar = 0.0 + returns;
#endif
    
    return(bstar);
}



// Self-employment Assistance Program -- cost
double bstarcostfun(double returns, int tgridfoc) {
    double bstarcost, btheta;
    
#if OPT_policy_exp == 2
    btheta = min((1.0-taxrate)*rhostar*wstar*prod[0]*parascale_DRI,(1.0-taxrate)*rhostar*wstar*prod[tgridfoc]);
#endif
#if OPT_policy_exp != 2
    btheta = (1.0-taxrate)*rhostar*wstar*prod[tgridfoc];
#endif
    
    if(returns <= 0.0) {bstarcost = btheta;}
    if(returns > 0.0 & returns < (btheta)/DRIpar) {bstarcost = btheta - DRIpar*returns;}
    if(returns >= (btheta)/DRIpar){bstarcost = 0.0;}
    
#if OPT_policy_exp == 3
    bstarcost = (1.0-taxrate)*rhostar*wstar*prod[tgridfoc];
#endif
#if OPT_policy_exp == 5
    bstarcost = 0.0;
#endif
    
    return (bstarcost);
    
}


// Self-employment Assistance Program -- probability losing it
double muendoSEA(double returns, int tgridfoc) {
    return (min(mupar,max(0.0,mupar*(bstarcostfun(returns,tgridfoc))/((1.0-taxrate)*rhostar*wstar*prod[tgridfoc]))));
}


/***************************************************/
/** FUNCTION CONSUMPTION (OR MAX WEALTH) FUNCTION **/
/***************************************************/

// Consumption Worker //
double findconsW(const double xval, const int igridfoc, const int tgridfoc, const int ygridfoc){
    return ((1.0+rstar)*grid[igridfoc] + yprod[ygridfoc]*prod[tgridfoc]*wstar*(1.0 - taxrate) - lumpsumSEA + 0.00000001 - xval);
}

// Consumption Unemployed Long-run //
double findconsU(const double xval, const int igridfoc, const int tgridfoc, const int Utype){
    double  UI;
    if(Utype == 0){UI = prod[tgridfoc]*wstar*(1.0 - taxrate)*rhostar - lumpsumSEA;}
    if(Utype == 1){UI = 0.0;}
    return ((1.0+rstar)*grid[igridfoc] + UI + minpar + 0.00000001 - xval);   // assume that they don't pay any taxes.
}

// USE TO compute the focB and focR val //
double findconsRB(const double xval, const double coh){return (coh + 0.000001 - xval);}

// USE to compute focE -- Entrepreneur who is insolvent kept a fraction.
double findCOH_B(const int igridfoc, const int zgridfoc, const int tgridfoc, double invest, const double Rphi, const int Jinsured){
    
    double profit_bef_exempt, coh, reste_invest, profit_kept;
    
    //if(invest < -0.0000001){printf("mistake invest %f",invest);getchar();} //invest = max(0.0,invest);
    //invest = max(0.0,invest);
    
    profit_kept = min(profitE(statez[zgridfoc],tgridfoc,invest),0);
    coh         = max(((1 - court_cost)*invest + profit_kept - frac_repay*(invest-grid[igridfoc])),0.0);
    
#if SEA == 1
    if(Jinsured == 1){coh = coh + bstarfun(0.0, tgridfoc);}
#endif
    
    coh = coh - lumpsumSEA; //+ minpar;//+ minpar*prod[tgridfoc];
    if(coh < 0.0){coh = 0.0;}
    return coh;
}


// USE to compute focE -- consumption non-excluded entrepreneur who repays at the end.
double profitRfun(const int igridfoc, const int zgridfoc, const int tgridfoc, const double invest, const double Rphi){
    double indicatorax;
    if(invest >= grid[igridfoc]){indicatorax = 1.0;}else{indicatorax = 0.0;}
    return(profitE(statez[zgridfoc],tgridfoc,invest) - (Rphi-1.0)*(invest - grid[igridfoc])*indicatorax);      // CAUTION : 1 - Rphi = rstar + wedgerate
}
double findCOH_R(const int igridfoc, const int zgridfoc, const int tgridfoc, double invest, const double Rphi, const int Jinsured){
    double indicatorax;
    if(invest >= grid[igridfoc]){indicatorax = 1.0;}else{indicatorax = 0.0;}
    //if(invest < -0.0000001){printf("mistake invest %f",invest);getchar();} //invest = max(0.0,invest);
    //invest = max(0.0,invest);
    double tot_returns = profitRfun(igridfoc,zgridfoc,tgridfoc,invest,Rphi);
    
#if SEA == 1
    if(Jinsured == 1){tot_returns = bstarfun(tot_returns,tgridfoc);}
#endif
    
    double coh = tot_returns + (rstar)*(grid[igridfoc]-invest)*(1.0-indicatorax) + grid[igridfoc] - lumpsumSEA;
    if(coh < 0.0){coh = 0.0;}
    return coh;
}


// USE to compute focE --  consumption non-excluded entrepreneur who repays at the end.
double profitEhatfun(const int igridfoc, const int zgridfoc, const int tgridfoc, const double invest){
    return(profitE(statez[zgridfoc],tgridfoc,invest));
}
double findCOH_Ehat(const int igridfoc, const int zgridfoc, const int tgridfoc, double invest, const int Jinsured){
    
    //if(invest < -0.0000001){printf("mistake invest %f",invest);getchar();} //invest = max(0.0,invest);
    //invest = max(0.0,invest);
    double tot_returns = profitEhatfun(igridfoc,zgridfoc,tgridfoc,invest);
    
#if SEA == 1
    if(Jinsured == 1){tot_returns = bstarfun(tot_returns, tgridfoc);}
#endif
    
    double coh = tot_returns + (rstar)*(grid[igridfoc]-invest) + grid[igridfoc] - lumpsumSEA;
    if(coh < 0.0){coh = 0.0;}
    return coh;
}


// STRUCTURE FOR OPTIMIZATION ---

// optimal saving
struct paramsgoldenRB{
    int igridCOHST, tgridST, zgridST, JinsuredST;
    double searchWoutST, C_val;
    double *expE_WWE_ST, *expE_UnE_ST, *expE_WWNE_JNE_ST, *expE_UnNE_JNE_ST, *expE_WWE_JE_ST, *expE_UnE_JE_ST;
#if OPT_compute_C_flow == 1
    double *C_expE_WWE_ST, *C_expE_UnE_ST, *C_expE_WWNE_JNE_ST, *C_expE_UnNE_JNE_ST, *C_expE_WWE_JE_ST, *C_expE_UnE_JE_ST;
#endif
#if SEA == 1
    double *expE_UiE_ST, *expE_WWNE_JNEi_ST, *expE_UiNE_JNEi_ST, *expE_WWE_JEi_ST, *expE_UiE_JEi_ST;
#if OPT_compute_C_flow == 1
    double *C_expE_UiE_ST, *C_expE_WWNE_JNEi_ST, *C_expE_UiNE_JNEi_ST, *C_expE_WWE_JEi_ST, *C_expE_UiE_JEi_ST;
#endif
#endif
};

// optimal saving
struct paramsgolden{
    int tgridST, igridST, zgridST, JinsuredST, iterST, methodST, ygridST, Utype, excludedST;
    double *valueBnextST, *valueRnextST,  *saveBnextST, *saveRnextST, *searchBnextST, *searchRnextST, *valueEhatnextST,  *saveEhatnextST, *searchEhatnextST;
    double searchJoutST, searchWoutST, C_val;
    double *JNE_save, *JNE_search, *JNE_default, *JNE_val, rate_endo, *JE_save, *JE_val, *JE_search;
    double *expU_UiNE_JNE_ST,*expU_UnNE_JNE_ST,*expU_UiE_JE_ST,*expU_UnE_JE_ST,*expU_UiNE_ST,*expU_UnNE_ST,*expU_UiE_ST,*expU_UnE_ST,*expU_WWE_ST,*expU_WWNE_ST,*expU_WWNE_JNE_ST,*expU_WWE_JE_ST;
    double *expW_WWE_JE_ST,*expW_WWNE_JNE_ST,*expW_WWE_ST,*expW_WWNE_ST;
#if OPT_compute_C_flow == 1
    double *C_expU_UiNE_JNE_ST,*C_expU_UnNE_JNE_ST,*C_expU_UiE_JE_ST,*C_expU_UnE_JE_ST,*C_expU_UiNE_ST,*C_expU_UnNE_ST,*C_expU_UiE_ST,*C_expU_UnE_ST,*C_expU_WWE_ST,*C_expU_WWNE_ST,*C_expU_WWNE_JNE_ST,*C_expU_WWE_JE_ST;
    double *C_expW_WWE_JE_ST,*C_expW_WWNE_JNE_ST,*C_expW_WWE_ST,*C_expW_WWNE_ST,*C_EhatnextST,*C_BnextST, *C_RnextST;
#endif
#if SEA == 1
    double *expU_UiE_JEi_ST, *expU_UiNE_JNEi_ST,*expU_WWNE_JNEi_ST,*expU_WWE_JEi_ST;
    double *expW_WWE_JEi_ST, *expW_WWNE_JNEi_ST;
    double *valueBinextST, *valueRinextST, *saveBinextST, *saveRinextST, *searchBinextST, *searchRinextST, *valueEhatinextST, *saveEhatinextST, *searchEhatinextST;
#if OPT_compute_C_flow == 1
    double *C_expU_UiE_JEi_ST, *C_expU_UiNE_JNEi_ST,*C_expU_WWNE_JNEi_ST,*C_expU_WWE_JEi_ST;
    double *C_expW_WWE_JEi_ST, *C_expW_WWNE_JNEi_ST;
    double *C_BinextST, *C_RinextST, *C_EhatinextST;
#endif
#endif
};


//zbrent to look for the optimal search level (only s_w or s_e when entrepreneur or worker)
struct searchEorWstruct{
    double diffValueNextfoc, Psw;
    int tgridparXX;
};

//zbrent to look for the optimal search level (when UNEMPLOYED only)
struct searchfocWwhenE{
    int tgridparXX, typeUU, excludedU, igridparXX;
    double tempWWNEval, tempUnNEval, tempUiNEval, tempWWEval, tempUnEval, tempUiEval, tempWWNE_JNEval, tempWWE_JEval, tempUnE_JEval, tempUiE_JEval, tempUnNE_JNEval, tempUiNE_JNEval, searchJoutST;
#if SEA == 1
    double tempWWNE_JNEival, tempUiNE_JNEival, tempUiE_JEival, tempWWE_JEival;
#endif
};



// main files
#include "saveforPMmethod.cpp"
#include "values/SEARCHING_FUN.cpp"
#include "computesearch.cpp"
#include "VFI.cpp"
#include "SIMUL_10y.cpp"
#include "SIMULATION.cpp"
#if indexPM == 1
#include "PMmethod.cpp"
#endif
#include "equilibrium.cpp"



int main(int argc, char* argv[]){
    
    
    /**********  SET OBSERVED MOMENT **********/
    // global
    obsmoments[0] = 2.65;      // Capital - output ratio (/10)
    obsmoments[1] = 0.57;      // Bankruptcy rate (frac of entrepreneur)
    
    // labor market
    obsmoments[2] = 5.9;       // Entrepreneur exit rate (*100)
    obsmoments[3] = 8.8;       // Fraction of Entrepreneurs (*100)
    obsmoments[4] = 5.0;       // Unemployment rate (*100)
    obsmoments[5] = 20.0;      // Fraction of entrepreneurs previously unemployment (*100)
    
    // wealth income distribution
    obsmoments[6] = 7.0;       // Ratio median net worth E/W
    obsmoments[7] = 1;         // Fraction with negative income (*10)
    
    // transition btw occupations
    obsmoments[8]   = 1.08;     // invrt-U-shape
    obsmoments[9]   = 0.85;     // invrt-U-shape
    obsmoments[10]  = 1.08;     // invrt-U-shape
    
    
    // files for outcome
    snprintf(distfileJ, sizeof(distfileJ), "OUTPUT/distribution/distJ_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    snprintf(distfileJi, sizeof(distfileJi), "OUTPUT/distribution/distJi_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    snprintf(distnecessity, sizeof(distnecessity), "OUTPUT/distribution/distnecessity_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    snprintf(distnecessityW, sizeof(distnecessityW), "OUTPUT/distribution/distnecessityW_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    snprintf(distfilewealth, sizeof(distnecessity), "OUTPUT/distribution/distfilewealth_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    
    snprintf(pricepathfile, sizeof(pricepathfile), "OUTPUT/project_method/pricepath_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    snprintf(CEV_J, sizeof(CEV_J), "OUTPUT/project_method/CEVJ_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    snprintf(CEV_W_file, sizeof(CEV_W_file), "OUTPUT/project_method/CEVW_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    snprintf(CEV_U_file, sizeof(CEV_U_file), "OUTPUT/project_method/CEVU_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    snprintf(VEILTfile, sizeof(VEILTfile), "OUTPUT/project_method/VEILTfile_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    snprintf(CEV_TOT, sizeof(CEV_TOT), "OUTPUT/project_method/CEV_TOT_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    snprintf(CEV, sizeof(CEV), "OUTPUT/project_method/CEV_POLICY=%d_OPT_labor_demand=%d_Teq=%d_Weq=%d_Req=%d_rhostar=%f_pLTpar=%f.txt", OPT_policy_exp, OPT_labor_demand, Teq, Weq, Req, rhostar, pLTpar);
    
    
    /*** WEIGHTING MATRIX ***/
    for(int l=0.0; l<nbmoments; l++){covmat[l][l] = 1.0;}
    
    
    /******* IF CONTROL RANDOM SEARCH IS OFF *******/
    //    #if CRS == 0
    double *PARAMETERS;
    PARAMETERS = (double *) calloc((nbpara), sizeof(double));
    
    PARAMETERS[0] = betapar;
    PARAMETERS[1] = kappaE;
    PARAMETERS[2] = kappaW;
    PARAMETERS[3] = psiparW;
    PARAMETERS[4] = nupar;
    PARAMETERS[5] = court_cost;
    PARAMETERS[6] = g1;
    PARAMETERS[7] = g2;
    PARAMETERS[8] = g3;
    PARAMETERS[9] = sigz;
    PARAMETERS[10] = pz;
    
    
    // compute equilibrium
#if (OPT_vary_UI_benefit + OPT_vary_UI_dur + OPT_vary_UI_both) == 0
    ComputeEquilibrium(PARAMETERS);
#endif
    
    // loop over the UI system.
#if OPT_vary_UI_benefit == 1
    double value_UI_repl[3] = {0.1,0.3,0.7};
    for(int pp = 0; pp < 3; pp++){
        pLTpar  = 0.5;
        rhostar = value_UI_repl[pp];
        ComputeEquilibrium(PARAMETERS);
    }
#endif
    
    // loop over the UI system.
#if OPT_vary_UI_dur == 1
    double value_UI_dur[3] = {0.9,0.15,0.075};
    for(int pp = 0; pp < 3; pp++){
        rhostar = 0.4;
        pLTpar  = value_UI_dur[pp];
        mupar   = value_UI_dur[pp];
        ComputeEquilibrium(PARAMETERS);
    }
#endif
    
    // loop over the UI system.
#if OPT_vary_UI_both == 1
    double value_UI_repl_2[4] = {0.4,0.4,0.6,0.8};
    double value_UI_dur_2[4]  = {0.25,0.132,0.5,0.5};
    for(int pp = 0; pp < 4; pp++){
        rhostar = value_UI_repl_2[pp];
        pLTpar  = value_UI_dur_2[pp];
        mupar   = value_UI_dur_2[pp];
        ComputeEquilibrium(PARAMETERS);
    }
#endif
    
    
    printf("PROGRAM FINISHED"); getchar();
    //    #endif
    
    
    //
    //    /******* IF CONTROL RANDOM SEARCH IS ON *******/
    //    #if CRS == 1
    //    int countsobol;
    //
    //    FILE *soboloutfile;
    //
    //    double *Best_PARA, *param_max, *param_min, ini_val2[nbsobol][nbpara];
    //    Best_PARA = (double *) calloc((nbpara), sizeof(double));
    //
    //    // PARAMETERS MIN MAX ///
    //    param_max = (double *) calloc((nbpara), sizeof(double));
    //    param_min = (double *) calloc((nbpara), sizeof(double));
    //
    //
    //    param_min[0] = 0.9737;      param_max[0] = 0.9743;     // BETA
    //    param_min[1] = 0.24;        param_max[1] = 0.28;       // KappaE
    //    param_min[2] = 0.84;        param_max[2] = 0.865;      // KappaW
    //    param_min[3] = 2.2;         param_max[3] = 2.6;        // Psipar
    //    param_min[4] = 0.789;       param_max[4] = 0.798;      // nupar
    //    param_min[5] = 0.0237;      param_max[5] = 0.0253;     // bank_cost
    //    param_min[6] = 0.0679;      param_max[6] = 0.0684;     // g1
    //    param_min[7] = 0.0774;      param_max[7] = 0.0780;     // g2
    //    param_min[8] = 0.1019;      param_max[8] = 0.1026;     // g3
    //    param_min[9] = 0.42;        param_max[9] = 0.45;       // sigz
    //    param_min[10] = 0.865;      param_max[10] = 0.885;     // pz
    //
    //
    //
    //    /******** INITIALIZE SOBOL SEQUENCE ********/
    //    #if SOBOLCREATE == 1
    //        // generate sobol
    //        float *sobolinitial;
    //        int isobol[nbsobol];
    //        sobolinitial = i4_sobol_generate(nbsobol, nbpara, 1100);
    //
    //        // SELECT ONLY SOBOL WHERE the MAPPING G() is INCREASING //
    //        countsobol=0;
    //        for(int i = 0; i < nbsobol; i++) {
    //            if(param_min[6]*(1.0 - sobolinitial[nbpara*i + 6]) + param_max[6]*sobolinitial[nbpara*i + 6] <= param_min[7]*(1.0 - sobolinitial[nbpara*i + 7]) + param_max[7]*sobolinitial[nbpara*i + 7] && param_min[7]*(1.0 - sobolinitial[nbpara*i + 7]) + param_max[7]*sobolinitial[nbpara*i + 7] <= param_min[8]*(1.0 - sobolinitial[nbpara*i + 8]) + param_max[8]*sobolinitial[nbpara*i + 8]) {isobol[i] = 100; countsobol++;}else{isobol[i] = 0;}
    //        //    printf("%f %f %f %d %d", param_min[6]*(1.0 - sobolinitial[nbpara*i + 6]) + param_max[6]*sobolinitial[nbpara*i + 6], param_min[7]*(1.0 - sobolinitial[nbpara*i + 7]) + param_max[7]*sobolinitial[nbpara*i + 7], param_min[8]*(1.0 - sobolinitial[nbpara*i + 8]) + param_max[8]*sobolinitial[nbpara*i + 8], countsobol, isobol[i]);getchar();
    //        }
    //
    //        //SAVE Sobol
    //        char buffer [80];
    //        int numbersobol = 0;
    //
    //        sprintf(buffer,"INPUT/SMM_para/SOBOL_%d.out",numbersobol);
    //        strncpy(SOBOL, buffer, sizeof(SOBOL) - 1);
    //        SOBOL[sizeof(SOBOL) - 1] = 0;
    //
    //        soboloutfile=fopen(SOBOL, "w");
    //        setbuf ( soboloutfile , NULL );
    //        for(int i = 0; i < nbsobol; i++) {
    //            for(int j = 0; j < nbpara; j++) {
    //                if(isobol[i] >= 100) {fprintf(soboloutfile,"%20.15f\n", sobolinitial[i*nbpara + j]);}
    //            }
    //        }
    //        fclose(soboloutfile);
    //    #endif
    //
    //
    //    #if SOBOLCREATE == 0
    //        // OR load sobol
    //        double *sobolinitial;
    //        sobolinitial = (double *) calloc((nbsobol*nbpara), sizeof(double));
    //        readinput(sobolinitial, nbsobol*nbpara, "INPUT/SMM_para/SOBOL_0.out");
    //        #endif
    //
    //        // b) Compute the corresponding points from the sequence
    //        for(int i = 0; i < nbsobol; i++) {
    //            for(int j = 0; j < nbpara; j++) {
    //                ini_val2[i][j] = param_min[j] + sobolinitial[nbpara*i + j] * (param_max[j] - param_min[j]);
    //            }
    //        }
    //
    //        for(int j = 0; j < nbpara; j++){
    //            printf("%f \t", ini_val2[0][j]);
    //        }
    //        printf("\n");
    //
    //        // WITHOUT OMP //
    //        myCRS(param_min, param_max, ini_val2, Best_PARA, ComputeEquilibrium);
    //
    //        // WITH OMP //
    //        #if OMP == 2
    //        myCRSomp(param_min, param_max, ini_val2, Best_PARA, ComputeEquilibrium);
    //        #endif
    //
    //        for(int i = 0; i < nbpara; i++) {
    //           printf("%20.15f\t", Best_PARA[i]);
    //        }
    //
    //        printf("PROGRAM FINISHED"); getchar();
    //
    //    #endif
    
    return 0;
}
